#ifndef _GEODESICS_H_
#define _GEODESICS_H_
//This file implement the modified Euler Method to compute the geodesics after applying FMM on the manifold
//#include <boost/numeric/mtl/config.hpp>


#define DISTANCE_ERROR 0.002
//This enumeraton indicate whether the next step should go along the edge or across a face
enum POINT_TYPE
{
	FACE_POINT, EDGE_POINT, VERTEX_POINT, END_POINT
};

typedef mtl::dense2D<double>  MATRIX;
typedef mtl::Collection<MATRIX>::value_type SCALAR;
typedef mtl::dense_vector<SCALAR>            VECTOR;

bool LU_solver (const MATRIX & A, const VECTOR & b, VECTOR & solution);
void Huen (double coeff[], const Mesh::Point & x_i, Mesh::Point & x_j, double h);
void Euler (double coeff[], const Mesh::Point & x_i, Mesh::Point & x_j, double h);
void convert_coordinate (Mesh::Point tri_A[3], Mesh::Point tri_B[3], MATRIX & rotation, VECTOR & translation);
void crossproduct (const VECTOR & v1, const VECTOR & v2, VECTOR & product);
bool is_in_triangle(const VECTOR tri[3], const VECTOR & p, int & edge_idx);
//R (A + t) = B
void trans_A2B(VECTOR & A, VECTOR & B, MATRIX & r, VECTOR & t);
bool intersection(Mesh::Point point[4], Mesh::Point & p);

bool three_points_interpolation(Mesh::Point & start, Mesh::Point v[], double value[], double step, Mesh::Point & intsec, int & next_face);
bool six_points_interpolation (const Mesh::Point & start, Mesh::Point v[], double value[], double step, Mesh::Point & intsec, double & intsec_value, int & next_face);
void draw_geodesics (Mesh & mesh, const Mesh::Point & start, Mesh::VertexHandle start_point_handle,
					 Mesh::FaceHandle face_handle, Mesh::HalfedgeHandle start_point_on_heh, 
					 POINT_TYPE point_type,  const Mesh::Point & end, 
					 const std::vector<Mesh::FaceHandle> & end_faces, double step, std::list<Mesh::Point> & path);
void generate_merged_mesh_geodesics(Mesh & mesh, std::list<Mesh::Point> & path);
void generate_seperate_mesh_geodesics(Mesh & mesh, std::list<Mesh::Point> & path);
void generate_seperate_group_geodesics(Mesh & mesh,  std::list<std::list<Mesh::Point> > & path);
void generate_seperate_mesh_point_cloud(Mesh & mesh, std::list<Mesh::Point> & path);

#endif
